function [FBx, FBy] = VFcomposition(AFx, AFy, BFx, BFy, width, height)

[x,y] = ndgrid(0:height-1,0:width-1);
xd = x - (height/2); yd = y - (width/2);

AFTx = x + AFx;
AFTy = y + AFy;

BFTx = x + BFx;
BFTy = y + BFy;

FBx = zeros(height, width);
FBy = zeros(height, width);

xBas0 = floor(BFTx);
yBas0 = floor(BFTy);
xBas1 = xBas0 + 1;
yBas1 = yBas0 + 1;

xCom = BFTx - xBas0;
yCom = BFTy - yBas0;
perc0 = (1-xCom).*(1-yCom);
perc1 = (1-xCom).*yCom;
perc2 = xCom.*(1-yCom);
perc3 = xCom.*yCom;

%% limit indexes to boundaries
check_xBas0 = (xBas0<0) | (xBas0>(height - 1));
check_yBas0 = (yBas0<0) | (yBas0>(width - 1));
xBas0(check_xBas0) = 0;
yBas0(check_yBas0) = 0;
check_xBas1 = (xBas1<0) | (xBas1>(height - 1));
check_yBas1 = (yBas1<0) | (yBas1>(width - 1));
xBas1(check_xBas1) = 0;
yBas1(check_yBas1) = 0;

Bx0 = AFx(1 + xBas0 + yBas0*height);
Bx1 = AFx(1 + xBas0 + yBas1*height);
Bx2 = AFx(1 + xBas1 + yBas0*height);
Bx3 = AFx(1 + xBas1 + yBas1*height);

By0 = AFy(1 + xBas0 + yBas0*height);
By1 = AFy(1 + xBas0 + yBas1*height);
By2 = AFy(1 + xBas1 + yBas0*height);
By3 = AFy(1 + xBas1 + yBas1*height);

%% set outside pixel to zero
Bx0(check_xBas0 | check_yBas0) = 0;
Bx1(check_xBas0 | check_yBas1) = 0;
Bx2(check_xBas1 | check_yBas0) = 0;
Bx3(check_xBas1 | check_yBas1) = 0;

By0(check_xBas0 | check_yBas0) = 0;
By1(check_xBas0 | check_yBas1) = 0;
By2(check_xBas1 | check_yBas0) = 0;
By3(check_xBas1 | check_yBas1) = 0;

FBx = Bx0.*perc0 + Bx2.*perc2 + Bx1.*perc1 + Bx3.*perc3;
%FBx = reshape(FBx, width, height);

FBy = By0.*perc0 + By1.*perc1 + By2.*perc2 + By3.*perc3;
%FBy = reshape(FBy, width, height);

FBx = FBx + BFx;
FBy = FBy + BFy;

% FBx = FBx - x;
% FBy = FBy - y;

end